### Ethiopia sentinel site field indicator summaries & spatial predictions
# Supporting documentation at: ...
# M.Walsh (Apr. 2013)

### Set local working directory e.g.,
# setwd("/Users/markuswalsh/Documents/Ethiopia/Field data")

### Required R-libraries, to install use:
# install.packages("arm"), dependencies=TRUE)
require(arm)

### Load data
fdat <- read.table("ET_field_dat.csv", header=T, sep=",")

### Mixed model-based, site-level predictions
# Prevalence predictions of cultivated areas
cult.glmer <- glmer(Cult~1+(1|Site), family=binomial, data=fdat)
display(cult.glmer)

# Extract and plot site-level predictions and standard errors
cult.coef <- coef(cult.glmer)
cult.se <- se.coef(cult.glmer)
coefplot(cult.coef$Site[,1], cult.se$Site[,1], varnames=rownames(cult.coef$Site), xlim=c(-6,6), CI=2, cex.pts=1.1, main="")

# Site level prevalence predictions of visible signs of soil erosion (VSSE)
vsse.glmer <- glmer(cbind(VSSE,4-VSSE)~1+(1|Site), family=binomial, data=fdat)
vsse.coef <- coef(vsse.glmer)
vsse.se <- se.coef(vsse.glmer)
coefplot(vsse.coef$Site[,1], vsse.se$Site[,1], varnames=rownames(vsse.coef$Site), xlim=c(-6,6), CI=2, cex.pts=1.1, main="")

# Mean difference in VSSE prevalence, between cultivated and non-cultivated sites
vsse.glmer1 <- glmer(cbind(VSSE,4-VSSE)~Cult+(Cult|Site), family=binomial, data=fdat)
display(vsse.glmer1)
anova(vsse.glmer, vsse.glmer1) 

# Prevalence predictions of soil depth restrictions to 30 cm (DR30) and mean differences between cultivated and non-cultivated areas 
dr30.glmer <- glmer(cbind(DR30,4-DR30)~1+(1|Site), family=binomial, data=fdat)
dr30.coef <- coef(dr30.glmer)
dr30.se <- se.coef(dr30.glmer)
coefplot(dr30.coef$Site[,1], dr30.se$Site[,1], varna (mes=rownames(dr30.coef$Site), xlim=c(-6,6), CI=2, cex.pts=1.1, main="")
dr30.glmer1 <- glmer(cbind(DR30,4-DR30)~Cult+(Cult|Site), family=binomial, data=fdat)
anova(dr30.glmer, dr30.glmer1)

# Prevalence predictions and mean differences in woody vegetation cover score of >15% (WCS) between cultivated and non-cultivated areas
wcs.glmer <- glmer(cbind(WCS,4-WCS)~1+(1|Site), family=binomial, data=fdat)
wcs.coef <- coef(wcs.glmer)
wcs.se <- se.coef(wcs.glmer)
coefplot(wcs.coef$Site[,1], wcs.se$Site[,1], varnames=rownames(wcs.coef$Site), xlim=c(-6,6), CI=2, cex.pts=1.1, main="")
wcs.glmer1 <- glmer(cbind(WCS,4-WCS)~Cult+(Cult|Site), family=binomial, data=fdat)
display(wcs.glmer1)
anova(wcs.glmer, wcs.glmer1)

# Site-level predictions and mean differences in shrub densities (nShrubs) between cultivated & non-cultivated areas
nshrub.glmer <- glmer(nShrubs~Cult+(1|Site), family=poisson, data=fdat)
nshrub.coef <- coef(nshrub.glmer)
nshrub.se <- se.coef(nshrub.glmer)
coefplot(nshrub.coef$Site[,1], nshrub.se$Site[,1], varnames=rownames(nshrub.coef$Site), xlim=c(-6,6), CI=2, cex.pts=1.1, main="")
nshrub.glmer1 <- glmer(nShrubs~Cult+(Cult|Site), family=poisson, data=fdat)
anova(nshrub.glmer, nshrub.glmer1)

# Site-level predictions and mean differences in tree densities (nTrees) between cultivated & non-cultivated areas by site
ntree.glmer <- glmer(nTrees~1+(1|Site), family=poisson, data=fdat)
ntree.coef <- coef(ntree.glmer)
ntree.se <- se.coef(ntree.glmer)
coefplot(ntree.coef$Site[,1], ntree.se$Site[,1], varnames=rownames(ntree.coef$Site), xlim=c(-6,6), CI=2, cex.pts=1.1, main="")
ntree.glmer1 <- glmer(nTrees~Cult+(Cult|Site), family=poisson, data=fdat)
anova(ntree.glmer, ntree.glmer1)

